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ABSTRACT 


The terahertz (THz) region of the electromagnetic spectrum covers frequencies ranging 
from approximately 100 GHz to 10 THz. This region of the spectrum has not been fully 
utilized due to the lack of compact and efficient sources as well as detectors. The aim of 
the present research is to explore the use of thin metal films as high THz absorbing 
materials and determine their THz absorbing characteristics analytically as well as 
experimentally. The films are to be used in bi-material based suspended structure, which 
sense minute changes in temperature due to THz absorption via difference in thermal 
expansion coefficients. The amount of deformation can be measured by reflecting a light 
beam from the pixel structure to determine the THz power incident on it. During the 
initial phase, design and fabrication of efficient THz absorbing multi-layer stack was 
carried out. The stack consists of a dielectric Bragg reflector and a metal film. The Bragg 
reflector serves to maximize reflection in the visible spectral range (required by the 
readout scheme), while the metal serves to provide strong absorption of THz radiation. 
The absorption characteristics of the stack were simulated using COMSOL finite element 
modeling and theoretical analysis performed using Fresnel's equations and then compared 
with the corresponding experimental results. The measured absorption characteristics of 
the stack agree well with the analysis. In addition, the thickness of the metal film was 
optimized for maximizing the THz absorption of the stack. 
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I. 


INTRODUCTION 


Recently, a significant number of applications that use terahertz (THz) radiation 
have emerged through deeper study of the THz region of the electromagnetic spectrum 
(100 GHz to 10 THz). Its non-ionizing* nature and the property of being able to penetrate 
most of insulating materials (textiles, paper, etc.) while reflecting from metallic objects, 
makes it a good candidate for the leading role in applications involving detection of 
concealed objects and medical imaging. To name a few, in the area of homeland security 
THz imagers can be used to detect concealed metal objects or certain types of explosives 
due to their characteristics THz signatures. Moreover, in the field of medicine, THz 
imaging appears particularly suited for examining skin abnormalities such as skin cancer 
[1] and characterizing human tissues [2]. 

Real time imaging in THz spectral range has been achieved using an uncooled 
microbolometer infrared (IR) camera and a milliwatt-scale quantum cascade laser (QCL) 
[3]. The principle of radiation detection is based on temperature changes arising from 
absorption of THz by an absorbing layer within the pixel. In addition, if the material used 
to build the pixel exhibits a high temperature coefficient of resistance (TCR), a relative 
difference in its resistance will be governed by the following equation [4]; 

R = R„(l + aAT) 

where is the initial resistance, a is the temperature coefficient of resistance and AT 
is the temperature difference. 

As a result and by assuming that each pixel has been constructed in such a way 
that a current bias (F^) can be applied via a Wheatstone bridge, one can measure the 

changes in voltage (AF^^J across the bridge as [4]; 

cxV 

AF^yj=^^AT (assuming aATD 1) 

* THz photons do not carry enough energy to remove (i.e., ionize) an electron from an atom or 
molecule and, therefore, cannot inflict significant damage to human DNA, a fact that is opposed to X-ray 
imaging technique. 
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The disadvantage of this scheme that bias used for reading the resistance change 
generates heat (self heating) in the pixel which limits the integration time affecting the 
signal-to-noise ratio [5]. On the other hand, an alternative and more efficient read out 
technique is under development. This detection scheme, as shown in Figure 1, involves 
deformation of a bi-material based suspended structure to minute changes in temperature 
due to difference in thermal expansion coefficients of layers of material used to build the 
pixel [6]. The amount of deformation can be measured by reflecting light off from the 
pixel structure (see Figure 1) and observing the reflected light using a CCD camera to 
generate an image of THz power incident on individual pixels. Furthermore, the optical 
readout technique does not require any readout electronics. Consequently, self-heating 
effects get eliminated when at the same time longer integration times are now available 
and hence better signal to noise ratios [7]. 


THz source 



Blade ^ 
Wrapped 
In Plastic 


Figure 1. The readout process begins by recording a background image using the 
reflected visible light from the array placed into the vacuum cell, when the THz 
source is off. After the THz illumination, through the lens and the wrapped object, 
the respective pixels which absorb that radiation will deform. As a result, all of or 
portion of the light reflected from the pixel array will be blocked by the aperture 
and result in the different amount of light reaching the CCD. By subtracting the 
background image, one can obtain the image (shadow) of the non THz transparent 
object. 
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The pixel arehiteeture, as displayed in Figure 2, undergoes a deformation due to 
absorption of certain amount of radiation. This deformation is achieved by designing the 
detector pixel to have areas of bi-material layers. The bi-material layers usually consist of 
an insulator and a metal with different coefficients of thermal expansion. When the 
temperature increases, these layers expand at different rates causing the whole structure 
to deform as shown in Figure 3. It has been shown that the use of SiOi-Al combination 
can provide a significant large deflection when operates in IR wavelengths [7]. Having 
high sensitivity in mind, it is expected that this type of detector would be appropriate for 
detection of THz radiation. Note that the power density in THz spectral band in the 
thermal background is much smaller than that of the infrared. For achieving high 
sensitivity in THz frequencies, it is important to develop strong THz absorbing 
membrane to integrate onto the pixel. 


Absorber 

with Optical Bimaterial 



Isolation 

Regions 


Figure 2. Schematics of the Bi-material Pixel’s Optimized for THz Frequencies 
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Subdomain: Total dspi ace merit [|jm] Deformation: Displacement Max: 3.868 



Min: O 


Figure 3. Simulated Deformation of the Pixel in Figure 2 using COMSOL Finite 
Element Modeling Software 

In this thesis, THz absorption characteristics of the membrane (i.e., a multi-layer 
stack of thin films) will be simulated using COMSOL finite element modeling and by 
theoretical analysis done using Fresnel’s equations. The aim is to develop a model that 
can estimate the absorption properties of multilayer stacks in the THz region that can be 
used to design sub-wavelength thin films with the maximum percentage of absorption of 
THz radiation. The experimental absorption measurements carried out using Fourier 
transform infrared spectroscopy (FTIR) of a thin-film stack developed in collaboration 
with Oakridge National Laboratory will also be presented [8]. 
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II. THEORY AND MODELING 


A. THEORETICAL BACKGROUND 

Absorption of electromagnetic radiation in a thin film stack containing both 
dielectric and metal films can be calculated by determining transmission and reflection 
coefficients as described in [9]. The percentage of absorption is then estimated by 
subtracting the two coefficients from the unity. In the following, the procedure described 
in [9] to calculate transmission and reflection coefficients is outlined by initially 
assuming a simple, three-layer stack as illustrated in Figure 4. 
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Figure 4. Three Layer Structure used in Calculation of Transmission and Reflection 
Coefficients (from [9]) 


The bottom and the top layer are simply air while the middle layer can be any 
material. The total reflection and transmission coefficients including multiple reflections 
in the film can be expressed as [9 and 10]; 


and 


Y Y e 

'l,2 ^'2,3 ^ 

\ ^ Y r e 

^ ^ 'l,2'2,3 ^ 


^2 


h,2C,3 ^ 

l + fi2fi,3 


-m 


( 1 ) 

( 2 ) 
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where and 




are the Fresnel coefficients of reflection and transmission, 


respectively, between layers k and m, and (p^ is the phase change when propagating 
through layer k. Such a phase change is given by 

'In 

(3) 


where X is the free space wavelength, Uk is the index of refraction of layer k, is its 
thickness, and 6^. is the associated term of the Snell’s law. The index of refraction of a 

metal thin film can be described by Equation 4 [10]. The complex nature of the index of 
refraction reflects the absorptive nature of the metal layer. The absorption of the 
incoming radiation will depend on the conductivity (a) of the metal. 

n, = (1 - i)^cT/(4ns J) , (4) 

where / is the frequency of incident radiation. 


This approach can be generalized for a stack containing arbitrary number of films, 
as depicted in Figure 5. 


Dielectric medium 

n. 

es (oo) 

Gold or silver thin film 

/\ ^ 

e 4 (10 to 100 nm) 

Chromium thin film 

A 

e 3 (0 to 5 nm) 

Glass plate 

/\ >^2 

Q 2 (0 or 1 mm) 

Coupling Prism 


ei (oo) 


Figure 5. A Thin Film Stack Containing Arbitrary Number of Films (from [9]) 


The reflection and transmission coefficients can be obtained using the following 
recursive formulas [9 and 10]; 


and 


r. = ■ 




-2J<Pi 




-2Jn 


(5) 
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v t 

where ' and ' are the total eoefficients of refleetion and transmission respeetively from 
i-th layer and above. Note that by assuming that the very first and last layers to be air, the 
eounter ^ starts from the seeond layer from the top and goes until 2 (i.e., in Figure 5 
where there are 5 layers, ^ starts from 4). A computer program was developed (see 
appendices A and B) to calculate the total reflection (fraction of the incoming power 
v v ^ 

reflected) ^ ^ and transmission (fraction of the incoming power transmitted) 

COS 0 0 0 

- —t 2 t 2 * coefficients as a function of THz frequency, where and are the 

cos Oin 

angles of incidence and transmission, respectively. By subtracting the reflection and 
transmission coefficients from the unity, absorption in the stack can be obtained as 

COS 0 

A = \-r2r2* - f^t2t2^ . (7) 

cos din 


The multi-layer stacks used in our experimental measurements consisted of the 
total of eleven thin films deposited on a 500 pm thick silicon substrate, as shown in 
Figure 6. 


75 nm 
110 nm 


500 pm 


SitJ, 


S1O2 

SiN. 


S1O2 

AiN, 


S1O2 

SiN" 


'' , 


S1O2 

SiN, 


S1O2 

Cr 

Si 


Figure 6. Structure of the Multilayer Stacks used in the Analysis and Measurement 
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Two multilayer stacks consists of thin Cr film of thickness either 15 nm or 30 nm 
was first grown on Si substrate followed by five successive periods of SiO^ and Si-^N^ 
with thicknesses of 110 nm and 75 nm, respectively. The SiO^ and Si^N^ layers act as a 

Bragg reflector for optical readout of pixel deflection as well as impedance matching to 
reduce the reflection at the metal interface. The indices of refraction we used for the 


dielectric materials S'Oj and Si^N^ were 1.46 and 2.05 respectively. From equation (4), 

which gives the corresponding complex index of refraction for Cr, the direct dependence 
of the conductivity of the metal of interest is apparent. The conductivities of Cr layers 
have been measured using a 4-pomt probing method and found to be 6.75 x 10^ [S/m] 
and 7.2 x 10^ [S/m] for 15 nm and 30 nm layers, respectively. These values are consistent 
with the measurements by Laman et al. [11] which showed that the conductivity of thin 
films is much lower than that of a bulk material due to high surface scattering. The 
complex index of refraction for the Si substrate was calculated using the following 
equation [12]. 



( 8 ) 


Where is the relative dielectric constant of silicon and cr is the conductivity of 
the silicon substrate. For the substrate used for the structure, the conductivity is about 10 
[S/m]. 

Using equations (l)-(8), the transmission, reflection and absorption of the two 
stacks consists of Cr layers of 15 nm and 30 nm were calculated for frequencies from 2.5 
to 10 THz pm at an angle of incidence of 30 degrees. The results are displayed in Figures 
7-12. The oscillatory behavior of the coefficients is due to the Fabry-Perot effect 
associated with the 500 pm thick Si substrate. 
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Figure 7. Calculated Transmission Coefficient of the Multi-layer Stack with 15 nm 
of Cr Layer as a Function of TFlz Frequency 


Theoretical Transmission (30nm of Cr) 
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Figure 8. Calculated Transmission Coefficient of the Multi-layer Stack with 30 nm 
of Cr Layer as a Function of TFlz Frequency 
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Theoretical Reflection (15nm of Cr) 
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Figure 9. Calculated Reflection Coefficient of the Multi-layer Stack with 15 nm of 
Cr Layer as a Function of TFlz Frequency 


Theoretical Reflection (30nm of Cr) 
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Figure 10. Calculated Reflection Coefficient of the Multi-layer Stack with 30 nm of 
Cr Layer as a Function of TFlz Frequency 
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Theoretical Absorption (15nm of Cr) 
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Figure 11. Calculated Absorption of the Multi-layer Stack with 15 nm of Cr Layer as 

a Function of TFlz Frequency 


Theoretical Absorption (30nm of Cr) 
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Figure 12. Calculated Absorption of the Multi-layer Stack with 30 nm of Cr Layer as 
a Function of TFlz Frequency 


It can be seen from the Figures 11 and 12 that the calculated absorption is about 
30% for the stack with 15 nm Cr and 20% for the one with 30 nm Cr. This is primarily 
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due to the high resistivity of the thinner film which generates more Joule heating and 
hence higher absorption. Note that above calculations give total absorption in Cr layer 
and Si substrate. 

B. FINITE ELEMENT MODELING 

The absorption characteristics of the stacks were also simulated using COMSOL 
finite element modeling. The main advantage of COMSOL is that it provides a 
convenient way to determine absorption in individual layers in the stack which is 
important when several layers absorb the incident radiation. However, the analytical 
approach developed above provides a way to compare the validity of the simulation. 
Before modeling using COMSOL, it is important to set appropriate physics and boundary 
conditions in addition to the physical dimensions of the layer structure in the stack. A 
detailed guidance for these necessary steps is given below [12]. 

I. Drawing the Multi-layer Stack 

Before starting to draw the stack, an appropriate selection of the physics through 
the model navigator in COMSOL has to be done. The fn-plane waves/Hybrid-Mode 
Waves/Harmonic propagation in the RF module were used in the modeling. This specific 
physics offers the harmonic electromagnetic waves analysis. Then the stack has to be 
drawn using the thicknesses exactly as they are described in the section on theoretical 
background. The lateral dimensions do not play any role because they are regarded as 
infinite by using the periodic boundary conditions. In addition, three extra layers of air 
have to be added on top and two at bottom by extending the very top and bottom layers 
so as to be at least 1 pm thick as illustrated in Figure 13. These two layers act as 
absorbers for the radiation exists (i.e., transmitted and reflected) from the stack. 
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Figure 13. Schematics of COMSOL Finite Elements Model 

Moreover, two composite objects were created by combining (a) the two top 
layers and (b) the rest of the layers. Then a boundary pair at the shared boundary of the 
two composite objects was created. The reason for these actions will be obvious later. 

2. Subdomain Settings 

After the stack has been designed, including the additional layers, the assignment 
of the corresponding materials for each layer has to be done. For the dielectric films, we 
just set n = 1.46 for Si02 and n = 2.05 for Si^N^. The refractive indices of chrome and 

silicon substrate are given by equations (4) and (8), respectively. These can be 
incorporated using the global expressions in COMSOL. Moreover, the very bottom and 
top extended layers (1pm) should be set to PML’s with absorption along y direction. 


13 















3. Boundary Settings 

The boundaries are divided into three groups. The first group ineludes the very 
top, the very bottom and the boundary pair (between the two eomposite objeets). These 
boundaries should be set as perfeet magnetie conductors. For the second group, which 
includes the other inner boundaries, the default setting (continuity) should be left 
unchanged. The settings for the third group corresponding to boundaries along the two 
sides are shown in Figure 14. 



Figure 14. Settings of the Side Boundaries of the Stack 

Above, and k^, are the x and y components of the k-vector respectively and they can 
also defined k^. in global expressions as 


k^ = ksin(6*), 

(9) 

k = k cos(6’), 

(10) 
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In 


( 11 ) 


k = - 


lambdaO _rjweh 


where 0 is the angle of incidence of the plain electromagnetic wave on the stack. 

Now focusing on the boundary pair at the boundaries between the two composite 
objects that has been created, this pair has to be set to the port boundary condition and the 
input power to be 1 Watt. Moreover, an analytical expression for the plane wave is 
required as shown in Figure 15. 



Figure 15. Setting the Analytical Expression for the Plane Wave 


After having done that, an incident field can be applied at the source boundary, 
and will only propagate in one direction into the domain, where it may reflect back. In 
addition, any reflected wave will pass back through the pair boundary unimpeded. Note 
that the source boundary is colored cyan and the destination purple by default. The source 
has to be at the bottom of the destination. 
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4. Meshing the Stack and Solving 

In order to follow the Nyquist theorem [14], at least eight mesh elements need to 
be present for each wavelength length. Since we work with very long wavelengths 
compared to the film thicknesses, i.e., from 30 pm to 120 pm, that rule does not play a 
significant role to our modeling procedure. Moreover, even our model is not symmetric 
(since we include the 500 pm Si substrate) and a manual meshing through sweeping 
methods should be necessary, a simple automatic mesh refining will work. 

By solving the problem for an angle of incidence of 30 degrees and wavelengths 
from 30 pm to 120 pm, we then can reproduce the corresponding reflection and 
transmission power (nPoav_rfweh) by simply using the postprocessing function in 
COMSOL and integrating it over the inner boundaries of the top and bottom PMLs, 
respectively. 

5. Comparison of the Absorption Results 

After subtracting the reflected and transmitted powers of the stack from 1 W of 
incident power, the amount of absorbed power by the stack can be estimated. The square 
points in Figures 16 and 17 show the modeled absorption for 15 and 30 nm stacks, 
respectively as a function of THz frequency. The solid lines in Figures 16 and 17 show 
the calculated absorption indicating a good agreement between the two approaches. 
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Figure 16. Theoretical Versus Modeled Absorption of the Multi-layer Stack with 15 
nm of Cr as a Function of THz Frequency 
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Figure 17. Theoretical Versus Modeled Absorption of the Multi-layer Stack with 30 
nm of Cr as a Function of TFIz Frequency 
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III. MEASUREMENT OE THZ CHARACTERISTICS 


A. FUNDAMENTALS OF FTIR SPECTROSCOPY 

The main instrument for conducting all the experimental measurements required 
for the absorption of thin-film materials in the THz range, was the NEXUS 870 FTIR. 
The incident THz radiation on the sample undergoes (a) reflection, (b) absorption and (c) 
transmission. The instrument is capable of measuring the intensity of radiation passing 
through the sample or reflecting from it. The transmission and reflection coefficients can 
be determined by measuring the same quantities either without the same or using a 
perfect reflector (gold coated mirror), respectively for subtracting the background. This 
process allows for a quantitative analysis of the resulting absorption of any multi-layer 
stack of interest. 

1. The Electromagnetic Aspect of the FTIR 

The functionality of the FTIR is based on a Michelson interferometer as shown in 
Figure 18. The light from a globar source in the FTIR spectrometer is divided at the 
beamsplitter, then reflected at either the translating or fixed mirror, and finally 
recombined at the beamsplitter to proceed to the sample area and the detector [7]. 


fi.ved tninvr 


source 



translating 

mirror 


beamsplitter 


detector 


Figure 18. Typical Diagram of Michelson’s Interferometer used in FTIR 
Spectroscopy 
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If the incoming electromagnetic wave is represented by: 

E = cos{cot -iKKy) , (12) 

where E^ is the average electric field amplitude and k the corresponding modified 

wavenumber (' ]), the detected intensity, I(x), as a function of the translating 
A. 

mirror’s displacement, x, is given by [8] 

I(x) = I(/f)[l + cos(2;r/fx)]. (13) 

In addition, one of the advantages of the FTIR is that many different wavenumbers due to 
a broadband source can be analyzed simultaneously. Thus, by first integrating over the 
wavenumbers and then taking Fourier transform of it, one can obtain a formula for the 
detected intensity as a function of the modified wavenumber k as, 

I(A:) = j[/(x)-^^^^]cos(2;rAx)d'x . (14) 

The software package (Omnic) came with the NEXUS 870 FTIR spectrometer 
perform the analysis based on the formula in equation (14) for extracting the spectral 
information. 

2. Description of Some Basic FTIR Parameters 

Before conducting any experimental measurement, the adjustment of some basic 
parameters of FTIR is recommended [15]. To name a few, the gain, the resolution, the 
number of scans and the velocity of the moving mirror. The default setting for the gain is 
the autogain. The adjustment of parameters to enhance specific signals has to be done 
with caution to avoid the saturation of the detector. There may be cases when a manual 
adjustment of the gain is necessary. The resolution has to do with the how sharp peaks in 
the spectrum can be resolved. The sharper the peaks in the spectrum the higher resolution 
needed. The default setting for the resolution in wavenumbers is 4 [c«t *], which 
correspond to an average data spacing of 2 [cm ' ]. The number of scans in combination 
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with the velocity of the moving mirror can play a significant role for the improvement of 
the signal-to-noise ratio. The default settings for the number of scans and mirror’s 

velocity are 32 and 6.3 [ —], respectively. 

s 

B. DESCRIPTION OF THE MULTILAYER STACK 

The multi-layer stacks used in our experimental measurements consisted of eleven 
thin films in total deposited on a 500pm thick silicon substrate, as shown in Figure 6. The 
detailed description of the two multilayer stacks used was given in Chapter II. 

The experimental measurements of the transmission and reflection coefficients of 
the stacks were conducted by using the Fourier transform interferometer (FTIR). By 
subtracting the sum of coefficients of reflection and transmission from 1 for each 
wavelength, one can determine the experimental value for absorption. All measurements 
were carried out for frequencies from 2.5 -10 THz. 

C. MEASURMENT OF TRANSMISSION COEFFICIENT 

The experimental setup for measuring the transmission of a wafer is shown in 
Figure 19. Before placing the stack in the sample compartment of the FTIR, the 
background spectrum is collected. In other words, the light beam produced by the FTIR 
to simply pass through the air from the right to the left without the sample in place. After 
having collected the background, the sample was placed and repeated the measurement. 
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The steps of the experimental proeedure are shown sohematieally in Figure 20 (a) 
and (b). By then properly adjusting the rotational base, the ineident angle of light on the 
sample was set at 30 degrees (for reasons that will be diseussed under the refleetion 
measurement). In addition, since a relatively low transmission through the stacks was 
obtained by modeling (see Figures 10 and 11), it would be appropriate to reduce the 
speed of the moving mirror of the FTIR from the default value to increase the light 
throughput for attaining higher signal-to-noise ratios [15]. The rest of the parameters left 
unchanged to avoid saturation during the background measurement. 
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Figure 20. a) Experimental Setup to Colleet the Baekground; 
b) Experimental Setup to Colleet the Sample’s Transmission 

The measured transmission eoeffieients for the two staeks with 15 nm and 30 nm 
of Chrome layers as a funetion of TFIz frequeney in 2.5 to 10 TFIz range are shown in 
Eigures 21 and 22, respeetively. The measured speetra are in elose agreement with that 
of the theoretical shown Eigures 7 and 8, respectively. In order to determine the 
absorption in the stacks, it is necessary to measure the reflection coefficients which will 
be described next. 
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Figure 21. Experimental Transmission of the Multi-layer Staek of 15 nm of Cr as a 
Funetion of THz Frequeney 
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Figure 22. Experimental Transmission of the Multi-layer Staek of 30 nm of Cr as a 
Funetion of THz Frequeney 

D. MEASUREMENT OF REFLECTION COEFFICIENT 

The experimental setup used for measuring the refleetion of a wafer is depieted in 
Figure 23. The light beam from the FTIR was steered by plaeing two gold eoated mirrors 

on the base at a fixed angle of 60 degrees with respeet to their axis of symmetry. This 
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corresponds to a 30 degrees angle of incidence on the sample. The background 
measurement was dome using a gold coated silicon wafer which behaves as a perfect 
mirror with close to 100% reflection. At this point, it is obvious why the angle of 
incidence for the transmission measurement was 30 degrees. 



Figure 23. Experimental Setup for Measuring the Reflection of a Wafer 


Again, the successive steps of the experimental procedure are schematically 
illustrated in Figures 24 (a) and (b). Since the simulated reflection from the stacks (see 
Figures 9 and 10) was relatively high the default settings of the FTIR were used in the 
measurement. The background measurement was carried out using a gold coated Si 

wafer to subtract the background from the sample, as shown in Figure 24 (a). 
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a) 


Background Scan 




Figure 24. a) Experimental Setup to Collect the Background; 
b) Experimental Setup to Collect the Sample’s Reflection 

The measured reflection coefficients of the stacks with 15 nm and 30 nm of 
Chrome as a function of TFIz frequency are shown in Eigures 25 and 26. The measured 
spectra are again in good agreement with that of the theoretical as shown in Eigures 9 and 
10 . 
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Figure 25. Experimental Refleetion of the Multi-layer Staek of 15 nm of Cr as a 
Funetion of THz Frequency 
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Figure 26. Experimental Reflection of the Multi-layer Stack of 30 nm of Cr as a 
Function of THz Frequency 
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E. DETERMINATION OF ABSORPTION 


In order to determine the experimental absorption of the stacks, one can subtract 
the sum of reflection and transmission coefficients from 1 for each frequency. The 
respective experimental results for the absorption of the stacks of 15 nm and 30 nm of 
Chrome as a function of the corresponding THz frequency are shown in Figures 27 and 
28. Note that these spectra include the contribution from the substrate due to its 
conductivity. Thought the conductivity of the substrate is relatively small compared to 
that of the Cr layer, it thickness is very large. For applications involving 
microcantilevers, the film stack is used without the substrate and it is important to 
determine the amount of TFIz absorption in the stack. Flowever, it is possible to 
determine the contribution from each layer by using COMSOL since allows the 
estimation of power dissipation in different layers in the stack [16]. In the following 
chapter, a comparison between the measured and simulated data and how to optimize 
absorption will be presented. 
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Figure 27. Experimental Absorption of the Multi-layer Stack with 15 nm Cr layer as 
a Function of TFIz Frequency 
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Figure 28. Experimental Absorption of the Multi-layer Stack with 30 nm Cr layer as 
a Function of TFlz Frequency 
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IV. DISCUSSION AND ANALYSIS 


A. INITIAL COMPARISONS 

Until now we have obtained theoretical, simulation and experimental data for 
conducting a comparison to examine the accuracy of the developed COMSOL simulation 
model. Figures 29-34 show comparison of experimental and simulated transmission, 
reflection and absorption spectra of the two multi-layer stacks with 15 nm and 30 nm 
thick Cr layers. Note that the simulation data has been reproduced with data spacing of 

2 [cm ] for better accuracy. In addition, the scales have been adjusted properly so as the 
comparisons to be as much distinct as possible. 
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Figure 29. Experimental Versus Simulated Transmission Coefficient of the Multi¬ 
layer Stack with 15 nm thick Cr as a Function of TFIz Frequency 
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Figure 30. Experimental Versus Simulated Transmission Coeffieient of the Multi¬ 
layer Staek with 30 nm thiek Cr as a Funetion of TFIz Frequeney 
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Figure 31. Experimental Versus Simulated Refleetion Coeffieient of the Multi-layer 
Staek with 15 nm thiek Cr Fayer as a Funetion of THz Frequeney 
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Figure 32. Experimental Versus Simulated Refleetion Coeffieient of the Multi-layer 
Staek with 30 nm thiek Cr Layer as a Funetion of TFlz Frequeney 
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Figure 33. Experimental Versus Simulated Absorption of the Multi-layer Staek with 
15 nm Cr Layer as a Funetion of TFIz Frequeney 
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Figure 34. Experimental Versus Simulated Absorption of the Multi-layer Stack with 
30 nm Cr Layer as a Function of TFIz Frequency 


It can be seen from the Figures 29-34 that the measurements agrees reasonably 
well with the finite element simulations done using COMSOL. As mentioned before, the 
oscillatory behavior seen in simulated data as well as to some extent in the measurements 
is due to Fabry-Perot effect associated with the Si substrate. In addition, Tables 1 and 2 
show comparisons of measured and simulated average values of transmission, reflection 
and absorption over the 2.5 to 10 THz range for the two multi-layer stacks with 15 nm 
and 30 nm thick Cr layers. The close agreement between the measurement and simulation 
shows that modeling can be used for optimizing the thickness of the Cr layer for strong 
THz absorption. 
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Table 1. Comparison of Average Coeffreients for the Staek with 15 nm thiek Cr 
Layer in the Spectral Region between 2.5 and 10 THz 


STACK of 15nm 
ofCr 

Transmission 

Coefficient 

(%) 

Reflection Coefficient 
(%) 

Absorption 

(%) 

Experimental 
(Average Value) 

5.87 

60.06 

34.07 

Respective standard 
deviation (%) 

0.21 

3.6 

3.7 

Model 

(Average Value) 

5.98 

60.02 

34.00 

Respective standard 
deviation (%) 

0.45 

6.1 

6.5 


Table 2. Comparison of Average Coefficients for the Stack with 15 nm thick Cr 
Layer in the Spectral Region between 2.5 and 10 THz 


STACK of 
30nm 
ofCr 

Transmission 

Coefficient 

(%) 

Reflection 

Coefficient 

(%) 

Absorption 

(%) 

Experimental 
(Average Value) 

4.36 

71.20 

24.55 

Respective 
standard 
deviation (%) 

0.22 

2.2 

2.1 

Model 

(Average Value) 

3.70 

73.14 

23.16 

Respective 
standard 
deviation (%) 

0.99 

4.1 

5.1 
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B. ANALYSIS AND OPTIMIZATION 


According to Tables 1 and 2, one can notice that a higher reflection resulted when 
the fllm thickness is increased from 15 to 30 nm. This is primarily due to the increase in 
conductivity of the Cr when the thickness is increased. In addition, Tables 1 and 2 also 
show a significant decrease of average absorption when shifting from 15 nm to 30 nm 
thickness of Cr layer. Although the transmission decreases, the reflection increases more 
and as a result the absorption drops. By inspection, one can notice that the increase in 
reflection is much larger than the transmission decrease. This suggests that the rate of 
decrease of transmission is smaller than the corresponding rate of increase of reflectance. 
In other words, as we go to thicker layers of Cr, the average absorption should 
proportionally decrease. 

As it has already mentioned in previous sections, the conductivities of Cr layers 
were measured to be 6.75 x 105 [S/m] and 7.20 x 105 [S/m] for 15 nm and 30 nm layers, 
respectively. Assuming a linear dependence of conductivity with the Cr layer thickness 
simulations can be performed to determine the optimum layer thickness for developing a 
stack with strong THz absorption. Based on the conductivity data of the two films, each 
2 nm increase of Cr layer thickness should represent a 0.046 x 105 [S/m] increase of 
conductivity. 

Since the thin-film stacks employed in THz bi-material pixels will not have the 
thick Si substrate, it is useful to optimize the absorption of the stack without the substrate. 
In addition, the angle of incident of THz on the pixel is normal than the 30 degrees used 
in earlier simulations to compare with experiments. Figures 35^2 show simulated 
absorption for a set of Cr layer thickness in 1 to 30 nm range as a function of THz 
frequency. Note that the oscillatory behavior disappeared from the simulated spectra due 
to removal of the substrate from the model. 
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Figure 35. Simulated Absorption of the Multi-layer Stack (without Substrate) with 1 
nm thick Cr Layer and Conductivity of 6.428 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 36. Simulated Absorption of the Multi-layer Stack (without Substrate) with 5 
nm thick Cr Layer and Conductivity of 6.520 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 37. Simulated Absorption of the Multi-layer Stack (without Substrate) with 8 
nm thick Cr Layer and Conductivity of 6.589 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 38. Simulated Absorption of the Multi-layer Stack (without Substrate) with 10 
nm thick Cr Layer and Conductivity of 6.635 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 39. Simulated Absorption of the Multi-layer Staek (without Substrate) with 11 
nm thiek Cr Layer and Conduetivity of 6.658 x 10^ [S/m] as a Funetion of TFlz 
Frequeney 
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Figure 40. Simulated Absorption of the Multi-layer Stack (without Substrate) with 13 
nm thick Cr Layer and Conductivity of 6.704 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 41. Simulated Absorption of the Multi-layer Stack (without Substrate) with 23 
nm thick Cr Layer and Conductivity of 6.934 x 10^ [S/m] as a Function of THz 
Frequency 
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Figure 42. Simulated Absorption of the Multi-layer Stack (without Substrate) with 30 
nm thick Cr Layer and Conductivity of 7.095 x 10^ [S/m] as a Function of TFlz 
Frequency 
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Table 3 shows the averaged values of simulated absorption of the staek over 2.5 
to 10 THz range as a funetion of Cr layer thiekness along with the eorresponding 
eonduetivity values used for simulation. Figure 43 shows a plot of absorption as a 
funetion of Cr layer thickness in the 1 to 30 nm range. From the plot in Figure 43, it is 
obvious that the optimum absorption occurs somewhere between 6 nm to 10 nm thick Cr 
layer with absorption of about 50% of the incident THz power. 


Table 3. Simulated Absorption for the Stack without the Substrate. The Middle 
Column gives the Conductivities used in the Simulation 


Varying thickness of Cr 
[nm] 

Reduced Conductivity of Cr [S/m] 

Modeled Absorption 
(%) 

1 

6.428 X 10^ 

19.10 

5 

6.520 X 10^ 

47.46 

7 

6.566 X 10^ 

50.29 

8 

6.589 X 10^ 

50.64 

9 

6.612 X 10^ 

50.55 

10 

6.635 X 10^ 

50.16 

11 

6.658 X 10^ 

49.55 

13 

6.704 X 10^ 

47.94 

15 

6.750 X 10^ 

46.06 

23 

6.934 X 10^ 

38.41 

30 

7.095 X 10^ 

32.83 
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Optimization of Absorption 
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Figure 43. Simualted Absorption of the Multi-layer Staek (without Substrate) as a 
Funetion of Chrome Layer Thiekness 
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V. CONCLUSION 


In the field of terahertz imaging, the use of MEMS bi-material sensors allows the 
implementation of optieal readout proeess whieh ean result better signal-to-noise ratio. In 
addition, it ean avoid the integration of readout eleetronies resulted in a simplified 
fabrieation proeess. For improving the sensitivity of bi-material pixels for THz imaging, 
it is neeessary to develop high THz absorbing thin films. In this thesis, an aeeurate model 
for simulating the absorptive behavior of thin-film materials in the THz range was 
developed using finite element modeling. The aeeuraey of the model was eompared with 
ealeulation done based on the analytieal approaeh deseribed in [9 and 10]. The two 
approaehes gave nearly identieal absorption speetra whieh eonfirmed the validity of the 
modeling. In addition, two thin-film staeks eontaining 15 and 30 nm Cr layers were 
eharaeterized using FTIR speetroseopy in 2.5 to 10 THz range. It was found that the 
strength of the absorption deereased from 34 % to 25 % when the film thiekness was 
inereased from 15 to 30 nm. This is primarily due to inerease in eonduetivity of the thin 
film whieh makes thieker film more refleetive. The measured absorption speetra agreed 
well with that of the simulations indieating the aeeuraey of the model for deseribing the 
absorption in thin film s . 

The model was used for optimizing the thiekness of the Cr layer for aehieving 
strong THz absorption. It was found that a Cr layer of approximately 10 nm gives nearly 
50 % absorption of ineident THz power. Though the model was developed for Cr metal 
films, ean be generalized to any metal layer of known eonduetivity. By then taking 
aeeount the proper funetionality of the developed model, the mierofabrieation of 
optimized absorptive layers of sub-wavelength thieknesses needed for use in MEMS bi¬ 
material deteetors, is aehievable. 

As has previously been stated in this thesis, a signifieant improvement in THz 
imaging ean be aehieved through the respeetive improvement in the readout teehnique. 
Based on a reeent researeh by Thomas Maier [17], mierobolometers ean be modified by 
metallie resonant absorber elements leading to an enhaneed responsivity at seleetable 
wavelengths. In addition, the dissipative energy absorption of tailored metamaterials 
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allowed for engineering the response of eonventional bolometer mierobridges. The 
absorption peak position and strength eould be determined by the geometry of the 
metamaterial. For square-shaped metal/dieleetrie metamaterial strueture showed speetral 
resonanees at wavelengths between 4.8 pm and 7.0 pm in aeeordanee with numerieal 
simulations using again COMSOL finite element modeling [17] with 80% absorption. 
The overall results suggest that there eould be also a similar application into the THz 
regime and could also be an interesting objective of a future research. 
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APPENDIX A: MATLAB CODE DEVELOPED EOR 
CALCULATING THE TOTAL POWER REELECTION OE THE 

MULTI-LAYER STACK 


nl = l; 
n2=2.05; 
n3=l.46; 


%Refraction Index of air 
% Refraction Index of Silicon Nitride layer 
% Refraction Index of Silicon Oxide layer 


n4=2.05; n5=1.46; n6=2.05; n7=1.46; n8=2.05; n9=1.46; nl0=2. 

nll=l.46; 


nl2=(1-i).*4500.*sqrt(lambda); %Refr. Index of Cr layer, eqn. (4) 
nl3=3.42-i.*108.6.*lambda; %Refr. Index of Silicon substrate,eqn. (8) 
nl4=l; 


d2=0.075e-6; % Thickness of Silicon Nitride layer 

d3=0.11e-6; % Thickness of Silicon Oxide layer 

d4=0.075e-6; d5=0.11e-6; d6=0.075e-6; d7=0.11e-6; 
d8=0.075e-6; d9=0.11e-6; dl0=0.075e-6; dll=0.11e-6; 

dl2=0.015e-6; % Thickness of Chrome layer 

dl3=500e-6; % Thickness of Silicon Substrate 


thetal=pi/6; 

w=8300:200:33300; 

"6 

lambda=l./w; 


% Angle of Incidence % 

% wavenumber range in meters, i.e., 2.5-lOTHz 

% wavelength, as expressed in FTIR % 


% Snell's law for every boundary of the multi-layer stack % 

theta2=asin(nl*sin(thetal)/n2);theta3=asin(nl*sin(thetal)/n3) ; 
theta4=asin(nl*sin(thetal)/n4);theta5=asin(nl*sin(thetal)/n5) ; 
theta6=asin(nl*sin(thetal)/n6) ; theta7=asin(nl*sin(thetal)/n7) ; 
theta8=asin(nl*sin(thetal)/n8);theta9=asin(nl*sin(thetal)/n9) ; 
thetalO=asin(nl*sin(thetal)/nlO);thetall=asin(nl*sin(thetal)/nil) ; 
thetal2=asin(nl.*sin(thetal) ./nl2);thetal3=asin(nl.*sin(thetal) ./nl3) ; 
thetal4=asin(nl.*sin(thetal)./nl4); 
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o\o o\o o\o LO o\o o\o o\o o\o o\o o\o 


% Change of phase in every layer of the multi-layer stack % 

phi2=2*pi./lambda*n2*d2*cos(theta2);phi3=2*pi./lambda*n3*d3*cos(theta3) ; 
phi4=2*pi./lambda*n4*d4*cos(theta4);phi5=2*pi./lambda*n5*d5*cos(thetaS); 
phi6=2*pi./lambda*n6*d6*cos(thetaG);phi7=2*pi./lambda*n7*d7*cos(theta7) ; 
phi8=2*pi./lambda*n8*d8*cos(theta8);phi9=2*pi./lambda*n9*d9*cos(theta9) ; 
phil0=2*pi./Iambda*nl0*dl0*cos(thetalO); 
phill=2*pi./lambda*nll*dll*cos(thetall) ; 
phi12=2.*pi./lambda.*nl2.*dl2.*cos(thetal2) ; 
phi13=2.*pi./lambda.*nl3.*dl3.*cos(thetal3); 


% Fresnel Coefficients of Reflection for every boundary % 

rl20=(nl*cos(thetal)-n2*cos(theta2))/(nl*cos(thetal)+n2*cos(theta2)); 
r23=(n2*cos(theta2)-n3*cos(theta3))/(n2*cos(theta2)+n3*cos(theta3)) ; 
r34=(n3*cos(theta3)-n4*cos(thetal))/(n3*cos(theta3)+n4*cos(thetal)) ; 
r45=(n4*cos(thetal)-n5*cos(thetaS))/(n4*cos(thetal)+n5*cos(thetaS)) ; 
r56=(n5*cos(thetaS)-n6*cos(thetaG))/(n5*cos(thetaS)+nG*cos(thetaG)) ; 
rG7=(nG*cos(thetaG)-n7*cos(theta7))/(nG*cos(thetaG)+n7*cos(theta7)) ; 
r78=(n7*cos(theta7)-n8*cos(theta8))/(n7*cos(theta7)+n8*cos(theta8)) ; 
r89=(n8*cos(theta8)-n9*cos(theta9))/(n8*cos(theta8)+n9*cos(theta9)); 
r910=(n9*cos(theta9)-nlO*cos(thetalO))/ 

(n9*cos(theta9)+nlO*cos(thetalO)) ; 
rl011=(nl0*cos(thetalO)-nll*cos(thetall))/ 

(nlO*cos(thetalO)+nll*cos(thetall)) ; 
rlll2=(nil.*cos(thetall)-nl2.*cos(thetal2))./ 
(nll.*cos(thetall)+nl2.*cos(thetal2)) ; 
rl213=(nl2.*cos(thetal2)-nl3.*cos(thetal3))./ 

(nl2.*cos(thetal2)+nl3.*cos(thetal3)) ; 
rl314=(nl3.*cos(thetal3)-nl4.*cos(thetal4))./ 

(nl3.*cos(thetal3)+nl4.*cos(thetall) ) ; 


% Total amplitude reflection for every effective layer % 

rl3=(rl213+rl314.*exp(-21.*phil3))./(1+r1213.*rl314.*exp(-21.*phil3)); 

rl2=(rlll2+rl3.*exp(-21.*phil2)) ./(1+r1112.*rl3.*exp(-21.*phil2)) ; 

rll=(rl011+rl2.*exp(-21.*phill))./(1+r1011.*rl2.*exp(-21.*phill)); 

rl0=(r910+rll.*exp(-2i*phil0)) ./ (l+r910*rll.*exp(-2i*phil0)); 

r9=(r89+rl0.*exp(-2i*phi9)) ./ (l+r89*rl0.*exp(-2i*phi9)); 

r8=(r78+r9.*exp(-2i*phi8)) ./ (l+r78*r9.*exp(-2i*phi8)); 

r7=(rG7+r8.*exp(-2i*phi7)) ./ (l+rG7*r8.*exp(-2i*phi7)); 

rG=(rSG+r7.*exp(-2i*phiG))./(l+rSG*r7.*exp(-2i*phiG)); 

rS=(r4S+rG.*exp(-2i*phiS)) ./ (l+r4S*rG.*exp(-2i*phiS)); 

r4=(r34+rS.*exp(-2i*phi4)) ./ (l+r34*rS.*exp(-2i*phi4)); 

r3=(r23+r4.*exp(-2i*phi3)) ./ (l+r23*r4.*exp(-2i*phi3)); 

r2=(rl20+r3.*exp(-2i*phi2)) ./(l+rl20*r3.*exp(-2i*phi2)) ; 


% TOTAL POWER REFLECTION COEFFICIENT % 


R2=r2.*conj(r2) ; 
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APPENDIX B: MATLAB CODE DEVELOPED EOR 
CALCULATING THE TOTAL POWER TRANSMISSION OE THE 

MULTI-LAYER STACK 


nl = l; 
n2=2.05; 
n3=l.46; 


%Refraction Index of air 
% Refraction Index of Silicon Nitride layer 
% Refraction Index of Silicon Oxide layer 


n4=2.05; n5=1.46; n6=2.05; n7=1.46; n8=2.05; n9=1.46; nl0=2. 
nll=l.46; 

nl2=(1-i).*4500.*sqrt(lambda); %Refr. Index of Cr layer, eqn.(4) 
nl3=3.42-i.*108.6.*lambda; %Refr. Index of Silicon substrate, eqn . (8) 
nl4=l; 


d2=0.075e-6; 
d3=0.lle-6; 


% Thickness of Silicon Nitride layer 
% Thickness of Silicon Oxide layer 


d4=0.075e-6; d5=0.lle-6; d6=0.075e-6; d7=0.lle-6; 
d8=0.075e-6; d9=0.lle-6; dl0=0.075e-6; dll=0.lle-6; 


dl2=0.015e-6; 
dl3=500e-6; 


% Thickness of Chrome layer 
% Thickness of Silicon Substrate 


thetal=pi/6; 

w=8300:200:33300; 

"6 

lambda=l./w; 


% Angle of Incidence % 

% wavenumber range in meters, i.e., 2.5-lOTHz 

% wavelength, as expressed in FTIR % 


% Snell's law for every boundary of the multi-layer stack % 

theta2=asin(nl*sin(thetal)/n2);theta3=asin(nl*sin(thetal)/n3); 
theta4=asin(nl*sin(thetal)/n4);theta5=asin(nl*sin(thetal)/n5) ; 
theta6=asin(nl*sin(thetal)/n6);theta7=asin(nl*sin(thetal)/n7) ; 
theta8=asin(nl*sin(thetal)/n8);theta9=asin(nl*sin(thetal)/n9) ; 
thetalO=asin(nl*sin(thetal)/nlO) ; thetall=asin(nl*sin(thetal)/nil) ; 
thetal2=asin(nl.*sin(thetal)./nl2);thetal3=asin(nl.*sin(thetal)./nl3); 
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o\° o\° o\° o\° o\° o\° (_n o\o o\° o\° 


thetal4=asin(nl.*sin(thetal)./nl4); 


% Change of phase in every layer of the multi-layer stack % 

phi2=2*pi./lambda*n2*d2*cos(theta2);phi3=2*pi./lambda*n3*d3*cos(theta3); 
phi4=2*pi./lambda*n4*d4*cos(theta4);phi5=2*pi./lambda*n5*d5*cos(thetaS) ; 
phi6=2*pi./lambda*n6*d6*cos(thetaG);phi7=2*pi./lambda*n7*d7*cos(theta7) ; 
phi8=2*pi./lambda*n8*d8*cos(theta8);phi9=2*pi./lambda*n9*d9*cos(theta9) ; 
phil0=2*pi./Iambda*nl0*dl0*cos(thetalO) ; 
phill=2*pi./lambda*nll*dll*cos(thetal1) ; 
phi12=2.*pi./lambda.*nl2.*dl2.*cos(thetal2); 
phil3=2*pi./lambda.*nl3.*dl3.*cos(thetal3) ; 


% Fresnel Coefficients of Reflection for every boundary % 

rl20=(nl*cos(thetal)-n2*cos(theta2))/(nl*cos(thetal)+n2*cos(theta2)) ; 
r23=(n2*cos(theta2)-n3*cos(theta3))/(n2*cos(theta2)+n3*cos(theta3)); 
r34=(n3*cos(theta3)-n4*cos(thetal))/(n3*cos(theta3)+n4*cos(thetal)) ; 
r45=(n4*cos(thetal)-n5*cos(thetaS))/(n4*cos(thetal)+n5*cos(thetaS)) ; 
r56=(n5*cos(thetaS)-n6*cos(thetaG))/(n5*cos(thetaS)+nG*cos(thetaG)) ; 
rG7=(nG*cos(thetaG)-n7*cos(theta7))/(nG*cos(thetaG)+n7*cos(theta7)); 
r78=(n7*cos(theta7)-n8*cos(theta8))/(n7*cos(theta7)+n8*cos(theta8)) ; 
r89=(n8*cos(theta8)-n9*cos(theta9))/(n8*cos(theta8)+n9*cos(theta9)) ; 
r910=(n9*cos(theta9)-nlO*cos(thetalO) ) / 

(n9*cos(theta9)+nlO*cos(thetalO)) ; 
rl011=(nl0*cos(thetalO)-nll*cos(thetall) ) / 

(nlO*cos(thetalO)+nll*cos(thetall)) ; 
rlll2=(nil.*cos(thetal1)-nl2.*cos(thetal2) ) . / 
(nll.*cos(thetall)+nl2.*cos(thetal2)) ; 
rl213=(nl2.*cos(thetal2)-nl3.*cos(thetal3))./ 
(nl2.*cos(thetal2)+nl3.*cos(thetal3)) ; 
rl314=(nl3.*cos(thetal3)-nil.*cos(thetal!))./ 

(nl3.*cos(thetal3)+nl4.*cos(thetal!)) ; 


% Total amplitude reflection for every effective layer % 

rl3=(rl213+rl314.*exp(-21.*phil3))./(l+rl213.*rl314.*exp(-21.*phil3)); 

rl2=(rlll2+rl3.*exp(-21.*phil2))./(l+rlll2.*rl3.*exp(-21.*phil2)); 

rll=(rl011+rl2.*exp(-21.*phill)) ./(l+rl011.*rl2.*exp(-21.*phill)) ; 

rlO=(r910+rll.*exp(-2i*phil0))./(l+r910*rll.*exp(-2i*phil0)); 

r9=(r89+rl0.*exp(-2i*phi9)) ./ (l+r89*rl0.*exp(-2i*phi9)); 

r8=(r78+r9.*exp(-2i*phi8)) ./ (l+r78*r9.*exp(-2i*phi8)); 

r7=(rG7+r8.*exp(-2i*phi7)) ./ (l+rG7*r8.*exp (-2i*phi7)); 

rG=(rSG+r7.*exp(-2i*phiG)) ./ (l+rSG*r7.*exp(-2i*phiG)); 

rS=(r4S+rG.*exp(-2i*phiS))./(l+r4S*rG.*exp(-2i*phiS)); 

r4=(r31+rS.*exp(-2i*phi4))./(l+r34*rS.*exp(-2i*phil)); 

r3=(r23+r4.*exp(-2i*phi3)) ./ (l+r23*r4.*exp(-2i*phi3)); 

r2=(rl20+r3.*exp(-2i*phi2)) ./ (l+rl20*r3.*exp(-2i*phi2)); 
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% Fresnel Coefficients of Transmission for every boundary % 

tl20=(2*nl*cos(thetal))/(nl*cos(thetal)+n2*cos(theta2)) ; 
t23=(2*n2*cos(theta2))/(n2*cos(theta2)+n3*cos(theta3) ) ; 
t34=(2*n3*cos(theta3))/(n3*cos(theta3)+n4*cos(theta4)); 
t45=(2*n4*cos(theta4))/(n4*cos(theta4)+n5*cos(thetaS)); 
t56=(2*n5*cos(thetaS))/(n5*cos(thetaS)+n6*cos(thetaG) ) ; 
t67=(2*n6*cos(thetaG))/(nG*cos(thetaG)+n7*cos(theta7)); 
t78=(2*n7*cos(theta7))/(n7*cos(theta7)+n8*cos(theta8)); 
t89=(2*n8*cos(theta8))/(n8*cos(theta8)+n9*cos(theta9)) ; 
t910=(2*n9*cos(theta9))/(n9*cos(theta9)+nlO*cos(thetal0)); 
tl011=(2*nl0*cos(thetalO))/(nl0*cos(thetalO)+nll*cos(thetall)) ; 
tlll2=(2.*nll.*cos(thetall)) ./(nll.*cos(thetall)+nl2.*cos(thetal2)) ; 
tl213=(2.*nl2.*cos(thetal2))./(nl2.*cos(thetal2)+nl3.*cos(thetal3)); 
tl314=(2.*nl3.*cos(thetal3))./(nl3.*cos(thetal3)+nl4.*cos(thetal4)); 


% Total amplitude transmission for every effective layer % 

tl3=(tl213.*tl314.*exp(-1.*phil3))./(l+rl213.*rl314.*exp(-21.*phil3)); 
tl2=(tlll2.*tl3.*exp(-1.*phil2))./(l+rlll2.*rl3.*exp(-21.*phil2)); 
tll=(tl011*tl2.*exp(-1.*phill))./(l+rl011*rl2.*exp(-2i*phill)); 
tl0=(t910.*tll.*exp(-1.*phil0))./(l+r910.*rll.*exp(-21.*phil0)); 
t9=(t89.*tl0.*exp(-1.*phi9)) ./ (l+r89.*rl0.*exp(-21.*phi9)); 
t8=(t78.*t9.*exp(-1.*phi8)) ./(l+r78.*r9.*exp(-21.*phi8)) ; 
t7=(tG7.*t8.*exp(-1.*phi7)) ./ (l+rG7.*r8.*exp(-21.*phi7)); 
tG=(t5G.*t7.*exp(-1.*phiG)) ./ (l+r5G*r7.*exp(-2i.*phiG)); 
t5=(t45.*tG.*exp(-1.*phi5)) ./ (l+r45.*rG.*exp(-21.*phi5)); 
t4=(t34.*t5.*exp(-1.*phi4)) ./ (l+r34*r5.*exp(-2i.*phi4)); 
t3=(t23.*t4.*exp(-1.*phi3))./(l+r23.*r4.*exp(-21.*phi3)); 
t2=(tl20.*t3.*exp(-1.*phi2))./(l+rl20*r3.*exp(-2i.*phi2)); 


% TOTAL POWER TRANSMISSION COEFFICIENT % 

T2=(nl4.*cos(thetall) ./ (nl.*cos (thetal))) .*t2.*conj (t2); 
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